/*********************************************************************************
 *
 * M&M CFD Simulator
 * Authors: Mohamed Ibrahim & Mohamed Bamakhrama
 * Copyrights (c) 2007. All Rights Reserved.
 *
 ********************************************************************************/

#include <stdio.h>
#include <stdlib.h>

#include "helper.h"
#include "visual.h"
#include "init.h"
#include "uvp.h"
#include "sor.h"
#include "boundary_val.h"
#include "freeboundary.h"
#include "cfd_utils.h"



/**
 * The main operation reads the configuration file, initializes the scenario and
 * contains the main loop. So here are the individual steps of the algorithm:
 *
 * - read the program configuration file using read_parameters()
 * - set up the matrices (arrays) needed using the matrix() command
 * - create the initial setup init_uvp(), init_flag(), output_uvp()
 * - perform the main loop
 * - trailer: destroy memory allocated and do some statistics
 *
 * The layout of the grid is decribed by the first figure below, the enumeration
 * of the whole grid is given by the second figure. All the unknowns corresond
 * to a two dimensional degree of freedom layout, so they are not stored in
 * arrays, but in a matrix.
 *
 * @image html grid.jpg
 *
 * @image html whole-grid.jpg
 *
 * Within the main loop the following big steps are done (for some of the
 * operations a definition is defined already within uvp.h):
 *
 * - calculate_dt() Determine the maximal time step size.
 * - boundaryvalues() Set the boundary values for the next time step.
 * - calculate_fg() Determine the values of F and G (diffusion and confection).
 *   This is the right hand side of the pressure equation and used later on for
 *   the time step transition.
 * - calculate_rs()
 * - Iterate the pressure poisson equation until the residual becomes smaller
 *   than eps or the maximal number of iterations is performed. Within the
 *   iteration loop the operation sor() is used.
 * - calculate_uv() Calculate the velocity at the next time step.
 * @author Mohamed Ibrahim and Mohamed Bamakhrama
 */
int main(int argc, char** args)
{

	/* The data structures needed during the simulation */
	struct geometry_data scenario_g;
	struct time_data scenario_t;
	struct pressure_iteration scenario_press;
	struct problem_parameters scenario_prob;
	struct boundary scenario_boundary;
	struct flow_vis scenario_vis;
	int response = 0;
	int n = 0;
	FILE *fp = NULL;
	int i, j;
	double dt_euler = 0.0;	/* used to solve Euler scheme in move_particles */

	nCount = 0;


	/* the initial value of res should be chosen carefully so that it enters the loop */
	scenario_press.res = 10.0;
	printf("*************************************************************************************************\n\n");
	printf("\t\t\t\tWelcome to M&M CFD Simulator v1.0\n\n");
	printf("*************************************************************************************************\n\n");
	printf("M&M: Starting up...\n");

	/* Check for the number of arguments supplied to the program through the command line */
	if (argc != 3) {
		/* print the usage and exit */
		printf("M&M CFD Simulator\n");
		printf("USAGE: sim [configuration_file] [output_file]\n");
		exit(-1);
	}
	else{

		printf("M&M: Initializing...\n");
		/* read the setup data from the configuration file */
		response = read_parameters( args[1],
					    &scenario_prob.Re,
					    &scenario_prob.UI,
					    &scenario_prob.VI,
					    &scenario_prob.PI,
					    &scenario_prob.GX,
					    &scenario_prob.GY,
					    &scenario_t.t_end,
					    &scenario_g.xlength,
					    &scenario_g.ylength,
					    &scenario_t.dt,
					    &scenario_g.dx,
					    &scenario_g.dy,
					    &scenario_g.imax,
					    &scenario_g.jmax,
					    &scenario_press.alpha,
					    &scenario_press.omg,
					    &scenario_t.tau,
					    &scenario_press.itermax,
					    &scenario_press.eps,
					    &scenario_t.dt_value,
					    &scenario_boundary.wl,
					    &scenario_boundary.wr,
					    &scenario_boundary.wt,
					    &scenario_boundary.wb,
					    &scenario_boundary.dp,
					    scenario_boundary.problem,
					    &scenario_vis.N,
					    &scenario_vis.x1,
					    &scenario_vis.y1,
					    &scenario_vis.x2,
					    &scenario_vis.y2,
					    &scenario_vis.dt_path,
					    &scenario_vis.dt_streak,
					    &scenario_vis.dt_insert,
					    &scenario_boundary.driving_velocity,
					    &scenario_boundary.ppc
					  );

		/* initialize the output file */
/*
		write_fieldfile( args[2],
				 scenario_g.imax,
				 scenario_g.jmax,
                 		 scenario_g.xlength,
				 scenario_g.ylength,
				 3
				);
*/
		scenario_t.t = 0.0;
		dt_euler = scenario_t.dt;

		/* allocate the needed memory */
		U = matrix ( 0 , scenario_g.imax + 1 , 0 , scenario_g.jmax + 1 );
		V = matrix ( 0 , scenario_g.imax + 1 , 0 , scenario_g.jmax + 1 );
		P = matrix ( 0 , scenario_g.imax + 1 , 0 , scenario_g.jmax + 1 );
		F = matrix ( 0 , scenario_g.imax + 1 , 0 , scenario_g.jmax + 1 );
		G = matrix ( 0 , scenario_g.imax + 1 , 0 , scenario_g.jmax + 1 );
		RS = matrix ( 0 , scenario_g.imax + 1 , 0 , scenario_g.jmax + 1 );

		Flag = imatrix( 0 , scenario_g.imax + 1 , 0 , scenario_g.jmax + 1 );
		init_imatrix( Flag, 0, scenario_g.imax + 1, 0 , scenario_g.jmax + 1, 0);

		Particle_Lines = init_particle ( scenario_boundary.problem,
						  &scenario_boundary.N,
						  scenario_boundary.ppc,
						  scenario_g.dx,
						  scenario_g.dy,
						  scenario_g.imax,
						  scenario_g.jmax,
						  U,
						  V,
						  &scenario_boundary.npart
						);
		/*fprintf(stderr,"%d\t%d\t\n",scenario_boundary.N, scenario_boundary.npart);*/
		if ( Particle_Lines != NULL ) {
			/*initial dt*/
			trace_fluid(	Particle_Lines,
					"fluid_particles",
					scenario_boundary.N,
					scenario_boundary.npart,
					scenario_t.t,
					scenario_g.xlength,
					scenario_g.ylength,
					dt_euler
				     );
		}

		init_uvp(	scenario_prob.UI,
				scenario_prob.VI,
				scenario_prob.PI,
				scenario_g.imax,
				scenario_g.jmax,
				U,
				V,
				P
			);
		/* TODO: Check if this has to be removed later on or not */

		init_flag( scenario_boundary.problem,
			   scenario_g.imax,
			   scenario_g.jmax,
			   Flag
			   );



		/*
		Streak_Lines = set_particle( scenario_vis.N,
					       scenario_vis.x1,
					       scenario_vis.y1,
					       scenario_vis.x2,
					       scenario_vis.y2
					       );
		Path_Lines = set_particle( scenario_vis.N,
					       scenario_vis.x1,
					       scenario_vis.y1,
					       scenario_vis.x2,
					       scenario_vis.y2
					       );
		*/
		/*
		mark_cells ( 	Flag,
				scenario_g.dx,
				scenario_g.dy,
				scenario_g.imax,
				scenario_g.jmax,
				scenario_boundary.N,
				Particle_Lines
				);
		*/

		/* Set the boundary values */
		boundaryvalues( scenario_g.imax,
				scenario_g.jmax,
				U,
				V,
				scenario_boundary.wl,
				scenario_boundary.wr,
				scenario_boundary.wt,
				scenario_boundary.wb,
				scenario_boundary.problem,
				Flag,
				scenario_boundary.driving_velocity
		    		);

		printf("M&M: Starting up the time-stepping loop...\n");

		while( scenario_t.t < scenario_t.t_end ) {

			calculate_dt( scenario_prob.Re,
				      scenario_t.tau,
				      &scenario_t.dt,
				      scenario_g.dx,
				      scenario_g.dy,
				      scenario_g.imax,
				      scenario_g.jmax,
				      U,
				      V
				      );

			mark_cells (	Flag,
					scenario_g.dx,
					scenario_g.dy,
					scenario_g.imax,
					scenario_g.jmax,
					scenario_boundary.N,
					Particle_Lines
					);

			/*if (scenario_t.t == 0 ) {*/
				fp = fopen("cells.txt", "w");
				for ( j = 0; j <= scenario_g.jmax; j++) {
					for ( i = 0; i <= scenario_g.imax; i++) {
						fprintf(fp, "%c", is_fluid(Flag[i][j]) ? 'F' : (is_obstacle(Flag[i][j]) ? 'O' : 'Y'));
					}
					fprintf(fp, "\n");
				}
				fclose(fp);
			/*}*/

			set_uvp_rand (	scenario_t.dt,
					scenario_g.dx,
					scenario_g.dy,
					scenario_g.imax,
					scenario_g.jmax,
					scenario_prob.Re,
					U,
					V,
					P,
					Flag,
					scenario_prob.GX,
					scenario_prob.GY
					);

			/* calculate F and G */
			calculate_fg( scenario_prob.Re,
				      scenario_prob.GX,
				      scenario_prob.GY,
				      scenario_press.alpha,
				      scenario_t.dt,
				      scenario_g.dx,
				      scenario_g.dy,
				      scenario_g.imax,
				      scenario_g.jmax,
				      U,
				      V,
				      F,
				      G,
				      Flag
				      );

			/* Calculate the right-hand side of the pressure eqaution */
			calculate_rs( scenario_t.dt,
				      scenario_g.dx,
				      scenario_g.dy,
				      scenario_g.imax,
				      scenario_g.jmax,
				      F,
				      G,
				      RS,
				      Flag
				      );


			/* Solve the system using w-SOR solver */
			scenario_press.it = 0;
			scenario_press.res = 2*scenario_press.eps; /* to gurantee that we enter the loop */
			while( scenario_press.it <= scenario_press.itermax && scenario_press.res > scenario_press.eps) {

				sor( scenario_press.omg,
				     scenario_g.dx,
				     scenario_g.dy,
				     scenario_g.imax,
				     scenario_g.jmax,
				     P,
				     RS,
				     &scenario_press.res,
				     Flag,
				     scenario_boundary.dp
				     );

				scenario_press.it += 1;
			}

			if ( scenario_press.it > scenario_press.itermax ) {
			/*	fprintf( stderr, "Warning: maximum number of iterations has been reached!\n");
			 */
			}


			/* update U and V with the new values */
			calculate_uv( scenario_t.dt,
				      scenario_g.dx,
				      scenario_g.dy,
				      scenario_g.imax,
				      scenario_g.jmax,
				      U,
				      V,
				      F,
				      G,
				      P,
				      Flag
				      );

			boundaryvalues( scenario_g.imax,
					scenario_g.jmax,
					U,
					V,
					scenario_boundary.wl,
					scenario_boundary.wr,
					scenario_boundary.wt,
					scenario_boundary.wb,
					scenario_boundary.problem,
					Flag,
					scenario_boundary.driving_velocity
			    		);

			set_uvp_rand (	scenario_t.dt,
					scenario_g.dx,
					scenario_g.dy,
					scenario_g.imax,
					scenario_g.jmax,
					scenario_prob.Re,
					U,
					V,
					P,
					Flag,
					scenario_prob.GX,
					scenario_prob.GY
					);
			dt_euler = fmin ( dt_euler, scenario_t.dt );
			move_particle (	Particle_Lines,
					scenario_boundary.N,
					dt_euler,
					scenario_g.dx,
					scenario_g.dy,
					scenario_g.imax,
					scenario_g.jmax,
					U,
					V,
					Flag
					);
			/*
			if ( Path_Lines != NULL ) {
				pathlines( Path_Lines,
					   scenario_vis.N,
					   "path_lines",
					   scenario_vis.dt_path,
					   scenario_t.t,
					   scenario_g.dx,
					   scenario_g.dy,
					   scenario_t.dt,
					   scenario_g.imax,
					   scenario_g.jmax,
					   U,
					   V,
					   Flag
					   );
			}
			*/
			if ( Particle_Lines != NULL ) {
				trace_fluid( Particle_Lines,
					     "fluid_particles",
					     scenario_boundary.N,
					     scenario_boundary.npart,
					     scenario_t.t,
					     scenario_g.xlength,
					     scenario_g.ylength,
					     dt_euler
					     );
			}

			/* Output the values used for visualization */
			if ( fabs (fmod(scenario_t.t, scenario_t.dt_value) - 0.0) < TOL ) {
				printf("M&M: Writing results to output files @ time step = %2.3f (iteration count = %d)\n",scenario_t.t,n);
				output_uvp( args[2],
					    scenario_g.xlength,
					    scenario_g.ylength,
					    scenario_g.imax,
					    scenario_g.jmax,
					    scenario_g.dx,
					    scenario_g.dy,
					    U,
					    V,
					    P,
					    Flag
					    );
			}
			scenario_t.t = scenario_t.t + dt_euler;/*scenario_t.dt;*/
			n += 1;
		}

	}

	/*
	printf("M&M: U[imax/2][7*max/8] = %5.5f\timax = %d\tjamx = %d\n",U[scenario_g.imax/2][7*scenario_g.jmax/8],scenario_g.imax,scenario_g.jmax);
	printf("M&M: Finalizing...\n");
	*/

	/* free the allocated memory */
	free_matrix( U,  0, scenario_g.imax + 1, 0, scenario_g.jmax + 1 );
	free_matrix( V,  0, scenario_g.imax + 1, 0, scenario_g.jmax + 1 );
	free_matrix( P,  0, scenario_g.imax + 1, 0, scenario_g.jmax + 1 );
	free_matrix( RS, 0, scenario_g.imax + 1, 0, scenario_g.jmax + 1 );
	free_matrix( F,  0, scenario_g.imax + 1, 0, scenario_g.jmax + 1 );
	free_matrix( G,  0, scenario_g.imax + 1, 0, scenario_g.jmax + 1 );

	free_imatrix( Flag,  0, scenario_g.imax + 1, 0, scenario_g.jmax + 1 );

	/* TODO: Free the memory used in the streaklines and pathlines */

	/* Yahooooo! we are done */
	printf("M&M: Simulation done...\n");
  	return 0;
}
